Gluon spectrum in the glasma from JIMWLK evolution 
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The JIMWLK equation with a "daughter dipole" running coupling is solved numerically, starting 
from an initial condition given by the McLerran-Venugopalan model. The resulting Wilson line 
configurations are then used to compute the spectrum of gluons comprising the glasma inital state 
of a high energy heavy ion collision. The development of a geometrical scaling region makes the 
spectrum of produced gluons harder. Thus the ratio of the mean gluon transverse momentum to the 
saturation scale grows with energy. Also the total gluon multiplicity increases with energy slightly 
faster than the saturation scale squared. 



PACS numbers: 24.85. +p,25.75.-q,12.38.Mh 
I. INTRODUCTION 

Particle production in the initial stage of a high en- 
ergy heavy collision, at RHIC or the LHC, is dominated 
by the small x gluonic degrees of freedom in the nuclear 
wavef unctions. When the dynamics of these gluons is 
dominated by a semihard saturation scale Q s it is possi- 
ble to understand the initial particle production in terms 
of weak coupling, first principles QCD. This can be done 
in the framework of the CGC effective theory (for re- 
views see e.g. pQ), where the calculation is organized in 
terms of the classical gluon field for the small x degrees 
of freedom, radiated from an effective color source repre- 
senting the larger x partons. This generic division only 
relies on the large energy, which, due to time dilation, 
makes the large x degrees of freedom evolve slowly, and 
on gluon saturation and weak coupling, which guarantees 
the validity of the classical field approximation. 

The energy (or x) dependence of the gluonic degrees of 
freedom is encoded in the JIMWLK 2] renormalization 
group equation, which is obtained by successively inte- 
grating out the quantum fluctuations around the classical 
field into the color source. The small x gluonic degrees of 
freedom are most conveniently described in terms of Wil- 
son lines in the classical field. The correlators of these 
Wilson lines can be probed experimentally by scatter- 
ing a dilute probe off the CGC, such as in DIS or pA 
collisions at forward rapidity. Many phenomenological 
applications of this framework use the BK [3] equation, 
which can be derived from JIMWLK in a mean field ap- 
proximation. 

A collision of two objects described in the CGC frame- 
work can also be understood in terms of classical gluon 
fields, the glasma fields 4 . They start off for r < 1/Q S 
as longitudinal chromoelectric and magnetic fields and 
evolve for r > l/Q s into modes that can be described as 
gluons with transverse momenta of order Q s . The basic 
structure of the glasma fields has been known for some 
time [5]. They can be obtained from a solution of the 
Classical Yang-Mills (CYM) equations of motion starting 
from an initial condition expressed in terms of the Wil- 
son lines mentioned above. The full equations have not 



been solved analytically, but there has been a significant 
amount of work to develop numerical solutions. So far, 
however, the numerical CYM [H[7] calculations have been 
performed using the Wilson lines from the MV model [8] , 
in stead of the full solution of the JIMWLK equation. 
The only calculations of gluon production to include the 
dynamics of high energy evolution have been done using 
a solution of the BK equation in a /cT-factorized approx- 
imation [5Hl2j. Since fey-factorization does not correctly 
give the initial gluon spectrum, or multiplicity, for the 
collision of two dense systems (the "AA" case) p!3Hl5] . 
these calculations are insufficient to give a full picture of 
the initial gluonic matter produced in a heavy ion colli- 
sion. 

The goal of this paper is to take the missing step 
of combining JIMWLK evolution with a full CYM cal- 
culation of gluon production in a heavy ion collision. 
The JIMWLK equation is solved numerically, using a 
("daughter dipole") running coupling constant. The re- 
sulting Wilson line configurations are then used as initial 
conditions in a numerical computation of the spectrum 
of gluons produced in a collision of two sheets of CGC. 
The point of view taken in this paper is that the MV 
model should be a reasonable initial condition for the 
evolution around RHIC energies. We shall show explic- 
itly how, starting from this initial condition, the effect of 
JIMWLK evolution is the creation of a geometrical scal- 
ing region for momenta kx > Q s , with a gluon spectrum 
that is harder than in the initial condition. This feature 
is carried over from the wavefunction to the spectrum 
of gluons in the glasma, making the glasma initial state 
more energetic at higher energies than a straightforward 
extrapolation of the MV model. 

The numerical method for solving the JIMWLK equa- 
tion is the one developed in [16; and the one for solving 
the CYM equations of motion that used e.g. in Refs 017]. 
We will thus only briefly describe them in Sec. [11], and re- 
fer the reader e.g. to Refs. [131 EH E] for a more exten- 
sive discussion. We will then characterize our results for 
the JIMWLK equation in Sec. [HI] and for the CYM cal- 
culation in Sec. |IV| before discussing phenomenological 
context and future directions in Sec. |V] 

The BK/ JIMWLK evolution is most conventionally 
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analyzed in terms of the fundamental representation sat- 
uration scale, which we denote by Q s . Gluon produc- 
tion, on the other hand, is expected to be dominated 
by the adjoint representation saturation scale which we 
denote Q s . These are related by a simple color factor 
Q 2 = [C A /C F }Q 2 = [2N C 2 /(N C 2 - 1)}QI The exact defi- 
nition of Q s used here is given in terms of the coordinate 
space Wilson line correlator in Eq. JTT1). 



II. 



SOLVING JIMWLK AND THE CYM 
EQUATIONS OF MOTION 



The JIMWLK equation describes the rapidity (or en- 
ergy) evolution of the probability distribution for Wilson 
lines. In the CGC framework large x degrees of free- 
dom are described by static color charges, which serve as 
sources for a classical color field. Most physical observ- 
ables can be expressed in terms of Wilson lines formed 
from the color field (in the covariant gauge, for a source 
moving in the positive z direction) 



U(xt) = -Pcxp J dx A+ ov (x T ,x )| 



(1) 



These Wilson lines are random SU(3) matrices from a 
probability distribution W y [U(xT)], which depends on 
the rapidity cutoff y that separates the large x color 
sources from the small x classical field. As the energy is 
increased, successive layers of quantum fluctuations have 
to be integrated into the probability distribution. This 
leads to the JIMWLK renormalization group equation 
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which describes the energy dependence of this probability 
distribution. This evolution equation can be written in 
a Langevin form for the rapidity evolution of the Wilson 
lines 



U y+dy (x. T ) = U y (x T )e 



ia°(x T ,!/)i° 



(3) 



where at each timestep the Wilson line is rotated in color 
space by 



a a (x T ,y) = - 
Ja s dy (x- z) 



ia s dy 
2tt 2 (xt — zt)'- 



Tr 
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with a random noise 

(r ?l a (x T , 2 ;)^(y T ,y')> = 5^S ab 5 2 ( XT - y T )S(y - v')- (5) 



. (4) 



Here U denotes the matrix in the adjoint representa- 
tion, and T a are the adjoint representation generators. 
This numerical procedure relies on the factorization of 
the JIMWLK kernel into a product of two terms, de- 
pending only on the coordinate pairs xt,zt and yt,zt 
(the "daughter dipoles"). 

Published numerical solutions of the JIMWLK equa- 
tion so far [THUS] have used a fixed coupling constant a s . 
This gives an evolution speed (increase of Q s with energy) 
that is too fast to be phenomenologically reasonable, so 
it is essential to use a running coupling constant here. 
There has been much discussion on the best running cou- 
pling prescription for BK/ JIMWLK evolution in the lit- 
erature 19 21j. The running coupling constant resums 
a subset of the NLO corrections to BK/ JIMWLK evo- 
lution, and different prescriptions correspond to resum- 
ming a different subset of them. Thus there is no unique 
"correct" way to set the scale of the coupling constant, al- 
though it has been argued [21] that the Balitsky prescrip- 
tion [50] minimizes the effect of other NLO corrections. 
The numerical Langevin method of solving the JIMWLK 
equation relies on factorizing the JIMWLK kernel into a 
product of two factors that only depend on the sizes of the 
"daughter" dipole x^ — z^ in Eq. Q. Thus implementing 
the preferred prescription of }20) would be difficult in a 
numerical solution of JIMWLK. We therefore use in this 
work the ad hoc "daughter dipole" prescription, where 
the magnitude of the coupling in Eq. Q only depends 
on xt — zt- One also has to regulate the Landau pole in 
the coupling constant, which we do in a smooth way at 
a scale /j,q by taking 



a 8 (r) 



12tt 



(33-2JV F )ln ( M 2 /A 2 ) 1/c + (r 2 A 2 /4) 



with N F = 3 and c = 0.2. 
pling is a — l2 ~" 



(6) 

The value of the frozen cou- 
The scale A in the cou- 



(33-2iV P )ln(Atg/A)- 

pling is parametrically of the order of Aqcd, but the 
exact value that should be used is scheme dependent. In 
the running coupling BK fit to DIS data [55] the scale is 
taken as a fit parameter and the data is found to prefer 
a smaller value A 2 = Aq CD /6.5, which we will assume 
here. 

The initial Wilson lines at y = are taken from 
the MV model, using the method discussed in more de- 
tail e.g. in Ref. [53]. After a given number of itera- 
tions of the Langevin equation ^ (we use a step size 
dy = 0.000l7r 2 /ao) one then obtains the configurations 
at a higher rapidity y. This is done separately for two 
independent configurations, corresponding to the two col- 
liding nuclei. In this work we only consider the symmetric 
situation of particle production at midrapidity, thus we 
evolve both nuclei starting from the same Q s o and for the 
same interval in y. 

Denoting the two nuclei as A and B one proceeds by 
constructing the light cone gauge fields corresponding to 



3 




FIG. 1: Wilson line correlator ( 10 1 in coordinate space; in 
lattice units (above) and as a function of the scaling variable 
rQ s (below). 



FIG. 2: Wilson line correlator ( |10| | (dipole cross section) in 
momentum space; as s function of krL (above) and of the 
scaling variable k/Q s (below). 



the Wilson lines as 



.4 



A.B 



-U\ B diU A .B- 

9 ' 



(7) 



The CYM equations of motion for the glasma fields are 
then solved as an initial value problem with initial con- 
ditions [5] 



A l \ 
A J1 \ 



=0 + 



=0+ 



— A / 
2 



(8) 
(9) 



At a given proper time of (r = 12/ Q s in our case) these 
fields are then Fourier-decomposed into k^-modes to de- 
termine the spectrum of the produced gluons. 



III. RESULTS FOR JIMWLK 

The most elementary observable to monitor during the 
evolution is the Wilson line correlator 

C(r = \x T - y T |) = -±r ( Tr C/t( XT )f/(y T )) . (10) 

The related dipole cross section 2 f d 2 b T (l — C(r)) ap- 
pears directly in the expression for the inclusive DIS cross 



section and is therefore of interest in itself. The correla- 
tor varies between the values 1 at r — and at large r, 
providing a natural way to define the saturation scale as 
the inverse of the correlation length of the Wilson lines 
as an intermediate scale between these two regimes. We 
adopt the definition, as in [24], of the fundamental rep- 
resentation Q s by the criterion 



C(r = V2/Q s ) = 



-1/2 



(11) 



Note that while for a Gaussian ("GBW"[5S]) correlator 
this definition is equivalent to the momentum space one 
used in Ref. [33], they need not give exactly the same 
values in the general case. 

In fcy-factorized calculations of gluon production in pA 
collisions one needs the unintcgratcd gluon distribution 
of the dense target. This is obtained from the the Fourier 
transform of Eq. ( 10 ) multiplied by k\ 



C(k 7 



d 2 r T e lkrrT C(r T ). 



(12) 



The typical behavior of the unintegrated distribution ( 12 ) 



is to start at zero for small kx and have maximum around 
kr ~ Qs- 

In the MV model the unintegrated gluon distribution 
behaves as C(kx) ~ \jk\ for large kr, which corresponds 
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FIG. 3: The evolution speed A = d In (x)/ d In 1/x as a 
function of the saturation scale Q s . The parameter values 
corresponding to the labels are detailed in Table [I] 



Configuration 


N± 


QsoL 




AL 


Base 


1024 


68 


15 


6 


Large Q s0 


1024 


136 


15 


6 


Large A 


1024 


68 


30 


12 


Small A 


1024 


68 


7.5 


3 


Large a 


512 


68 


15 


6 


Small Qo 


1024 


68 


30 


6 


Small L 


512 


34 


7.5 


3 


Small L, a 


1024 


34 


7.5 


3 



TABLE I: The values of the numerical parameters used in 
our simulation sets. For no /A — 2.5 the coupling is frozen at 
the value cto = 0.762 and for /Uo/A = 5 at Qo = 0.434. 

to the integrated gluon distribution xg(x,Q 2 ) behaving 
as - InQ 2 . The main effects of JIMWLK/BK evolution 
are the increase of the characteristic scale Qs with energy 
and and making the functional form of the unintcgratcd 
distribution less steep, ~ 1/kp . For fixed coupling the 
anomalous dimension is |26| 7 ~ 0.63 and for running 
coupling numerical solutions [27 s of the BK equation give 
7 ~ 0.85. Our present calculation is done on a linear (as 
opposed to logarithmic) lattice, and cannot go to very 
large momenta before lattice ultraviolet cutoff effects are 
felt in the spectrum. Therefore one cannot hope to get a 
very good numerical evaluation of the anomalous dimen- 
sion. The change in the behavior of the unintegrated 
gluon distribution is, however, clearly observable. The 
Wilson line correlators at different rapidities are shown 
in Fig. [T] in coordinate space and in Fig. [2] in momentum 
space, both in lattice units and a functions of the scaling 
variables rQ s and k/Q s . Both the increase of Q s and the 
development of a geometric scaling region are very well 
visible. 

The values used in the different sets of numerical com- 
putations in this paper are summarized in Table [I] The 
dependence of physical observables on the collision en- 
ergy and rapidity depends most of all on the speed of 




FIG. 4: Gluon spectrum at different energies, labeled by 
the rapidity interval of evolution starting from the MV initial 
condition at y — 0. The momentum is scaled by the saturation 
scale Q s corresponding to the rapidity in question. 

evolution, conventionally parametrized as 

A^ dlnQ ^. (13) 
dy 

At fixed coupling A ~ a s , and at running coupling the 
speed is expected to be proportional to a s (Q 2 ). Thus the 
evolution speed is controlled by the relation of the initial 
saturation scale to the QCD scale controlling the running 
of the coupling, Q s0 /A. In a physically realistic case we 
would like to start the evolution at a scale corresponding 
to midrapidity at RHIC energies, i.e. Q s o ~ 1.1 GeV [2"3"] . 
corresponding to Q s o/A ps 11, i.e. Q s o/A ps 7.6. Assum- 
ing the transverse area to be L 2 = 140 fm 2 this leads to 
QsqL w 68. These correspond to the "base" set of values 
in Table |T] and used in Figs, [I] [2]and|4j To test the effect 
of the initial scale on the speed of evolution we also show 
results with two different values of Q s o/A. The change in 
QsO I A is obtained by increasing Q s o by a factor of 2 and 
keeping all other parameters fixed. Alternatively A is de- 
creased or increased by a factor of 2, with everything else 
fixed. The value at which the coupling freezes without 
altering the dynamics at higher momentum scales can be 
altered by changing [i . The sensitivity of the calcula- 
tion to lattice effects has been tested by changing the 
lattice spacing a (i.e. the lattice UV cutoff ~ 1/a) and 
the physical volume L = Nj_a by a factor of two, with 
all the other dimensionful parameters fixed. 

The initial conditions in the MV model are con- 
structed [33] as a product of N y = 100 infinitesimal Wil- 
son lines, with the MV model color charge density param- 
eter g 2 /i adjusted to provide the desired saturation scale. 
The evolution speed for the different configurations is 
shown in Fig. [3j as a function of Q s /A. This confirms 
the expectation that starting with a lower Q s q/A results 
in a faster initial evolution. The phenomenologically pre- 
ferred evolution speed A < 0.3 is only reached with an 
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— — base 

• large Q sQ 

o large A 
□ large a 

a Q = 0.43 
o small L, a 

• small A 



y 



FIG. 5: Gluon liberation coefficient as a function of collision 
energy, parametrized by the rapidity interval of evolution from 
the initial scale. 



initial saturation scale that is higher than the above es- 
timates corresponding to RHIC energy. Note, however, 
that this is affected by the uncertainty concerning the 
correct value of A in this running coupling scheme. 



IV. RESULTS FOR THE GLUON SPECTRUM 

We then take the Wilson line configurations from the 
JIMWLK evolution and, as discussed in Sec. [IlJ use them 
as initial conditions for the solution of the CYM equa- 
tions of motion. The gluon spectrum resulting from the 
calculation is shown in Fig. [3] One can see that the 
gluon spectrum gets gradually harder as one moves from 
the initial condition into the geometric scaling regime. 
In the rapidity interval considered here the total multi- 
plicity and transverse energy of the gluons are still finite. 
This is to be contrasted with the case of fixed coupling 
where, for an unintegrated gluon distribution behaving 



as C(kx) 



■.-27 



the produced gluon spectrum would 



behave as dN/ d 2 px ~ Pj^"' ' ■ If the geometric scaling be- 
havior continued to arbitrary large kx this would result 
in an ultraviolet divergent transverse energy for 7 < 0.75, 
including the fixed coupling value ~ 0.63. In practice the 
geometrical scaling region does not extend up to infinite 
kx, but a very hard gluon spectrum could be difficult to 
reconcile with the observed transverse energy of the later 
stages of the quark gluon plasma. 

The dominant transverse momentum scale of the pro- 
duced gluon spectrum is expected to be the adjoint rep- 
resentation Q s . It is convenient to parametrize the gluon 
spectrum by the dimensionless "liberation coefficient" [251 
|2T)] c proportional to the total gluon multiplicity 



dN, 



init. g 



Cu 



d 2 x-r dy 2ir 2 c 



(14) 



^-base 
* large Q s0 
o large A 
□ large a 




and the mean transverse momentum of the produced glu- 
ons (pr) ■ For a more thorough discussion on relating 



FIG. 6: Mean gluonic transverse momentum as a function of 
collision energy.. 



these values to the measured charged particle multiplici- 
ties we refer the reader to e.g. [50] , 

The values of c and (pr) /Q s for different amounts of 
evolution are plotted in Figs. [5] and [6j We see that, as 
expected, the harder /c^-dependence in the initial wave- 
functions leads to a harder spectrum of produced gluons, 
as evidenced by the increase of (pr) /Q s from the initial 
condition. This increase does, however, seem to saturate 
after y ss 3. This indicates that the functional form starts 
to settle towards a new scaling characteristic of JIMWLK 
evolution, where the anomalous dimension 7 < 1 leads 
to a harder spectrum than in the MV model. The mean 
transverse momentum of the gluons in the glasma, scaled 
by Q s , increases from ~ 1 for the MV model initial con- 
dition to around ~ 1.5. 

A somewhat less expected feature is the increase of 
the scaled multiplicity c seen in Fig. [5] This is seen 
also in the shape of the scaled spectrum in Fig. |4j where 
new gluons are added for pr > 2Q S , but the shape for 
Pt < Qs changes very little. This observation could, as 
in Ref. [T3], be interpreted as a difference between "ini- 
tial" and "final" state interactions. At high transverse 
momentum pr > 2/Q s the produced gluon spectrum re- 
flects the unintegrated gluon distributions in the colliding 
projectiles. Thus fcr-factorization works in this regime, 
and the spectrum gets harder because of the develop- 
ment of the geometrical scaling window. The shape of 
the spectrum for pr < Q s , on the other hand, is a result 
of nonlinear interactions in the glasma stage, which ren- 
der the spectrum IR finite. These interactions are not 
captured in the fcr-factorized formalism. 

Both c and (pr) are smaller for the configurations 
where Q s a is large ("small a" and "large Q s "), which im- 
plies a strong dependence on the lattice ultraviolet cutoff 
l/a. This is to be contrasted with Fig. [3j where the evo- 
lution speed A showed no significant dependence on the 
UV cutoff. Testing this by further decreasing a with the 
same L would become prohibitively expensive for this 
study. We can, however, achieve smaller values of Q s a 
by decreasing the size of the system, i.e. Q S L. Making 
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FIG. 7: The Wilson line correlator after 5.2 units in rapidity. 
The points and the dashed line differ only by the value of 
the lattice spacing a used. The dotted and solid line likewise 
differ only by the value of a, and are obtained with L half of 
the first ones. For the parameter values used see Table [I] 

Q s qL smaller is eventually limited by the goal of staying 
in the strong field regime already for the initial rapidity. 

The effect of the lattice UV cutoff is demonstrated fur- 
ther in Figs. [7| and [8| They show the unintegrated gluon 
distribution (Fig. u¥ and the produced gluon spectrum 
(Fig. [8]) after 5.2 rapidity units of evolution for two dif- 
ferent values of Q s qL, with two different values of Q s qcl 
each, keeping the other parameters of the evolution the 
same. One can see that with the new phase space open- 
ing up at smaller a the unintegrated gluon distribution 
is mostly unaffected. The produced gluon spectrum, on 
the other hand, increases in the high px tail. Thus the 
produced gluon spectrum is more sensitive to the lattice 
regularization than the JIMWLK evolution itself. One 
consequence of this is that the leveling off of the increase 
in (pt)/Qs at high rapidity, seen in Fig. [6j could be in- 
fluenced by lattice cutoff effects. 

V. DISCUSSION 

We have, for the first time, used the solution of nonlin- 
ear the high energy evolution equations in a calculation of 
gluon production in the initial stage of a heavy ion colli- 
sion without resorting to a /cy-factorized approximation. 
This enables us to compute the gluon multiplicity and 
transverse energy density without ambiguities related to 
the normalization and without additional infrared cut- 
offs. It has been seen that the effect of JIMWLK evo- 
lution is to make the gluon spectrum harder, leading to 
a growth of the total multiplicity that is slightly faster 
than ~ and a gluon mean px that grows faster than 
Q s . The resulting gluon spectrum in the glasma, shown 
in Fig. [4] is the main result of this paper. 

The numerical calculation confirms that, as expected 
from studies of the BK equation, introducing a running 
coupling slows down the evolution to a speed more consis- 
tent with experimental observations. At fixed coupling, 




FIG. 8: Scaled gluon spectrum after 5.2 units in rapidity. 
The labels are as in Fig. [7] 

numerical JIMWLK evolution is known to be very sensi- 
tive to the lattice ultraviolet cutoff. With running cou- 
pling the speed of evolution becomes essentially indepen- 
dent of the UV cutoff. However, when the solution is 
used as an input for a calculation of the gluon spectrum 
of the initial glasma phase of heavy ion collisions, the 
lattice spacing dependence again becomes stronger. A 
full systematic continuum extrapolation is, however, left 
for future work. Effects of higher order corrections to 
JIMWLK/BK evolution could significantly modify the 
physics at kx ^ Qs- This would have a larger effect on 
the spectrum of gluons in the glasma than expected just 
from the speed of the evolution. 

The inclusion of NLO effects in our calculation has, 
by reasons of of numerical practicality, been limited to 
a "daughter dipole" prescription for the running cou- 
pling. Incorporating more of the NLO corrections to 
JIMWLK/BK evolution into the could have a much 
larger effect on the gluon spectrum in a nucleus-nucleus 
collisions than on, say, the total DIS cross section. An- 
other, separate but phenomenologically extremely topical 
issue that can be addressed in this same framework, are 
long range rapidity correlations in the glasma [TT1, 15T1 152"] . 
The calculation of the "ridge" correlation in the glasma 
proceeds in a similar fashion as the one performed in this 
paper, but is left for future work. 
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